1 Présentation du projet

Ce projet s’inscrit dans le cadre de la formation DSA, et évalue l’apprentissage de R réalisé auprès de Robin RYDER.

Le thème de ce projet est libre.

J’ai choisi d’étudier une série d’articles qui m’ont frappé sur la comparaison des approches classiques par modèles linéaires généralisés (GLM) et par machine learning (ci après ML) dans le cadre de la tarification en assurance non-vie et la manière dont elles pouvaient se compléter.

1.1 Les bonnes propriétés du GLM dans le cadre d’une modélisation de la fréquence de sinistres

Soit $N_{i}$ le nombre de sinistre de la police ${i}$, alors dans le modèle Poisson homogène on suppose $N_{i} \sim \mathscr{P}(\lambda \nu_{i})$ où $\nu_{i}$ est l’exposition au risque de la police $i$, $\lambda$ est la fréquence de sinistre annuelle (homogène sur le portefeuille)et $\mathscr{P}$ la loi de Poisson. On pose $Y_{i} = \frac{N_{i}}{\nu_{i}}$ la fréquence annuelle des sinistres pour la police $i$ et $ = (x_{1},x_{2},, x_{q}) $ les $q$ caractéristiques de la police $i$.

En tarification, on suppose dans les faits que $\lambda$ dépend de $\mathbf{x}$, de sorte que $\mathbb{P}[Y=k/\nu] = \mathbb{P}[N=k] = e^{-\lambda(\mathbf{x})\nu} \frac{(\lambda(\mathbf{x})\nu)^{k}}{k!}$ avec $ln(\lambda(\mathbf{x})) = \beta_{0} + \beta_{1}x_{1} + \beta_{2}x_{2} + \cdots + \beta_{q}x_{q} \overset{def}{=} \langle \mathbf{\beta}, \mathbf{x} \rangle$ en complétant $\mathbf{x}$ par $x_{0} = 1$.

Lorsqu’on utilise les approches classiques de GLM, on estime $\hat{\lambda}$ en utilisant l’estimateur du maximum de vraisemblance $\hat{\beta}$ de $\beta$.

Cet estimateur peut être obtenu soit en maximisant la log vraisemblance, soit de manière équivalente en minimisant la déviance de Poisson.

En complément, lorsqu’on utilise un GLM incluant un intercept et sa fonction de lien canonique, l’utilisatiuon d’un estimateur du maximum de vraisemblance apporte de bonnes propriétés pour son usage en tarification dont celle dite de propriété d’équilibre (Property 2.4 p33 de Data Analytics for Non-Life Insurance Pricing), qui précise que :

$\displaystyle\sum_{i=1}^{n} \nu_i\hat{\lambda}(\mathbf{x_{i}}) = \displaystyle\sum_{i=1}^{n} \nu_i exp\langle \mathbf{\hat{\beta}}, \mathbf{x_{i}} \rangle = \displaystyle\sum_{i=1}^{n} N_{i}$

De sorte que le nombre total de sinistres estimés est égale au nombre total de sinistres observés.

Cette propriété est même plus étendue dans certains modèles GLM, notamment lorsque le modèle poissonien avec lien log (exemple (b) p423 de A systematic relationship between minimum bias and gener- alized linear models) : dans ce cas l’équilibre existe au globaldu portefeuille, mais aussi pour chaque classe des variables catégorielles incluses dans le modèle.

Ces propriétés d’équilibres sont critiques pour un assureur dans un exercice de tarification, puisqu’elles assurent qu’une structure tarifaire basée sur cette approche lui permettra d’équilibrer (au moins théoriquement) ses encaissement avec les sinistres qu’il s’est engagé à indemniser.

Le GLM est donc interprétable, souple à l’usage et présente des propriétés d’équilibre global et local. Néanmoins son ajustement nécessite un temps humain important et le praticien doit s’en remettre à son expérience pour estimer les éventuelles interactions à inclure dans son analyse.

1.2 L’apport et les limites des modèles ML en tarification

Les modèles ML par leur grande flexibilité sont des candidats naturels pour remplacer ou à minima compléter les modèles GLM du fait de leur performance sur une grande variété de cas d’usages ces dernières années. Toutefois, ces modèles nécessitent qulques ajustements avant de pouvoir être utilisés dans des problèmes de tarification.

1.2.1 Ajustement #1 : définition de la fonction de perte

Une première adaptation des modèles classiques consiste donc à définir la déviance de poisson comme fonction de perte des algorithmes de ML.

C’est notamment l’approche préconisée par dans Data Analytics for Non-Life Insurance Pricing de Mario V. Wuthrich et Christoph Buser.

Cette approche a notamment conduit à l’article Boosting insights in insurance tariff plans with tree-based machine learning methods et la créatuon du package distRforest. NB : le package rpart inclut la déviance de poisson, mais ne permet de définir simplement une fonction de déviance appropriée pour la sévérité (lognormale ou gamma)

1.2.2 Limites de la modification de la fonction de perte

La définition d’une fonction de perte appropriée est devenue l’approche par défaut de l’utilisation des méthodes ML en tarification. Néanmoins, dès 2019 Wütrich remarqua que les sinistres estimés par les modèles ML ne correspondaient pas à la sinistralité observée sur le portefeuille, même sur le jeu d’entraînement

Il semble que l’utilisation directe de modèles de ML, mêmes adaptés pour minimiser la déviance, ne présentent pas les garanties d’équilibre des GLM.

Les modèles de ML permettent une meilleurs corrélation avec la réponse individuelle de chaque police, du fait de la plus grande flexibilité accordée à la forme de la contribution (éventuellement non linéaire) de chaque variable à la sortie. Toutefois, cette meilleure corrélation s’obtient au prix de l’abandon de la notion d’équilibre. En effet, aucune garantie n’est expressément introduite au global du portefeuille ce qui peut produire des divergences, potentiellement importantes, entre la somme des primes demandées et les sinistres que l’assureur s’est engagé à régler.

Cet effet est attribué par M Wütrich aux conditions d’arrêt précoce des algorithmes d’optimisation utilisés par les modèles ML.

1.3 Pistes de réconciliation des approches GLM et ML pour la tarification

Deux options s’offrent dès lors à l’actuaire :

essayer d’extraire des informations de la manière dont les modèles ML ont appris les efefts des variables sur la sortie pour compléter les GLM classiquement manipulés avec l’objectiof de combler en partie l’écart de performance (s’il est important) avec le modèle ML tout en conservant les propriiétés et l’interprétabilité naturelle des modèles GLM. C’est la voie suivie par l’article Peeking into the black box de Christian Lorentzen et Michael Mayer
corriger les approches ML en les combinant avec une approche GLM pour en récupérér les propriétés. Le gain de cette approche devrait être plus important que pour la méthode précédente tout en présentant les propriétés d’équilibre des modèles ML sans toutefois apporter l’interprétabilité de ces modèles. C’est la proposition de de Michel Denuit, Arthur Charpentier et Julien Trufin dans Autocalibration and Tweedie-dominance for Insurance Pricing with Machine Learning repartent du constat de Wütrich et proposent une solution de contournement sous la forme d’un ajustement en deux temps successifs :
1. Dans un premier temps, un modèle ML est ajusté pour déterminer une première estimation $\hat{\pi}$ des primes théoriques
2. Dans un deuxième temps, ces primes théoriques $\hat{\pi}$ sont utilisées comme variable entrante d’un modèle d’autocalibration local pour produire des primes corrigées $\widehat{\pi_{BC}}$. Cette méthode consiste à ajuster des modèles GLM locaux réduits à un intercept sur des ensembles de points “proches”. La notion de proximité ici se fait qui sont définis à partir de la proximité des primes théoriques estimées par ML $\hat{\pi}$ (l’étape de ML rassemble les observations qui se ressemblent, l’étape locale GLM lisse les écarts entre observations voisines pour retrouver les propriétés d’équilibre).

1.4 Objectif du projet

Dans ce rapport nous prendrons comme prétexte l’analyse du jeu de données tarifaire d’un portefeuille MRH qui nous avait été confié lors du Hackathon pour explorer ces deux approches. Ce faisant nous nous écarterons du cadre du Hackathon qui proposait de minimiser le RMSE. Par ailleurs, ne disposant pas du jeu de test, nous ne pourrons comparer ces méthodes avec le leaderbord produit à l’occasion de l’événement.

De ce fait, et sans perte de généralité, vis-à-vis du propos des articles, nous nous restreindrons dans la suite de ce document à l’analyse de la seule fréquence des sinistres.

Dans la section suivante nous présenterons le jeu de données avant d’en entreprendre la modélisation. La dernière section traitera des approches mises en places dans les articles.

2 Analyse

2.1 Mise en place de l’environnement

Le code suivant est un ensemble d’utilitaire pour naviguer dans les répertoires relatifs du projet et installer ou charger les packages nécessaires à l’analyse :

if (!require("here")){ 
  install.packages("here") 
  library("here")
}

set_here

LoadPackage <- function (load.lib=c("")) {

  install.lib<-load.lib[!load.lib %in% installed.packages()]
  for(lib in install.lib) install.packages(lib,dependencies=TRUE)
  sapply(load.lib,require,character=TRUE)
  
}


LoadPackage(c('ggplot2', 'dplyr', 'formattable', 'DT', 'tidyr', 'caret', 'plotly', 'xgboost', 'flashlight', 'MetricsWeighted', 'corrplot', 'lsr', 'h2o', 'wesanderson'))

2.2 code

plot_obs<-function(df,feature,use_year=T) {

    if (use_year) {
          dt <-df %>%
          group_by(df[[feature]], ANNEE, .drop=FALSE) %>%
          summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) ,  pct_EXPO = sum(EXPO) /  sum(df$EXPO)) %>%
          rename(feat = 'df[[feature]]') 
    } else {
          dt <-df %>%
          group_by(df[[feature]], .drop=FALSE) %>%
          summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) ,  pct_EXPO = sum(EXPO) /  sum(df$EXPO)) %>%
          rename(feat = 'df[[feature]]') 
    }
    
    
    
    fig <-plot_ly()
    fig <- fig %>% layout(title = paste0('Fréquence de sinistres par ', feature),
             xaxis = list(title = feature),
             yaxis = list(side = 'left', title = 'Exposition', showgrid = FALSE, zeroline = FALSE),
             yaxis2 = list(side = 'right', overlaying = "y", title = 'Fréquence', showgrid = FALSE, zeroline = FALSE))
    
    if (use_year) {
        for (year in unique(dt$ANNEE)) {
        
              dt_temp <- dt %>% filter(ANNEE==year)
              fig <- fig %>% add_trace( 
                             x= dt_temp$feat, y= dt_temp$EXPO, type = 'bar', name = paste0('Exposition_', year),
                             hovertemplate = 'Expo: %{y:.0f}<br>'
                          )
              fig<- fig %>% add_trace(
                             x= dt_temp$feat, y= dt_temp$Freq, type = 'scatter', mode = 'lines', yaxis = 'y2', name= paste0('Fréquence_', year),
                             hovertemplate = 'Freq: %{y:.2%}<br>'
                            )
        
        }
    } else {
              dt_temp <- dt 
              fig <- fig %>% add_trace( 
                             x= dt_temp$feat, y= dt_temp$EXPO, type = 'bar', name = 'Exposition',
                             hovertemplate = 'Expo: %{y:.0f}<br>'
                          )
              fig<- fig %>% add_trace(
                             x= dt_temp$feat, y= dt_temp$Freq, type = 'scatter', mode = 'lines', yaxis = 'y2', name= 'Fréquence',
                             hovertemplate = 'Freq: %{y:.2%}<br>'
                            )
      
    }
    
    return(fig) 
}

resume <- function(df, feature) {

    dt <-df %>%
      group_by(df[[feature]], .drop=FALSE) %>%
      summarise(Freq = sum(NB*EXPO)/sum(EXPO), nb_lignes = n(), EXPO=sum(EXPO) ,  pct_EXPO = sum(EXPO) /  sum(df$EXPO ))
    
    colnames(dt)[1] <- feature
    
    
    table<- formattable(dt,
                align = c("l", rep("r", NCOL(dt) - 1)),
                list(
                    Freq = percent,
                    nb_lignes= accounting,
                    EXPO = accounting,
                    pct_EXPO = percent
                )
    )
  
    fig1<- plot_obs(df, feature, use_year = F)
    
    fig2<- plot_obs(df, feature, use_year = T)
    
    return(list(dt=dt, table=table, fig=fig1, figyr=fig2))
}

3 Présentation du jeu de données

Le jeu de données est constitué de 4 fichiers :

list.files(path=here('data', 'raw'))

## [1] "expo_test.csv"  "expo_train.csv" "primes2019.csv" "sin_train.csv"

Le projet suit un portefeuille d’assurance MRH suivi sur plusieurs années de 2016 à 2018 inclus. L’objectif initial du Hackathon est d’estimer les primes 2019 sur un échantillon de contrats.

3.1 expo_train

Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :

expo_train = read.table(file = here('data', 'raw', 'expo_train.csv'), header=T, sep=',', dec='.', encoding = 'UTF-8', stringsAsFactors = T)

datatable(head(expo_train))

str(expo_train)

## 'data.frame':    155651 obs. of  19 variables:
##  $ X                  : int  384538 441079 119668 5124 130065 450830 151401 296526 71801 360445 ...
##  $ EXPO               : num  1 0.825 1 1 1 ...
##  $ FORMULE            : Factor w/ 4 levels "ALL_INCLUDE",..: 4 2 3 3 3 2 3 2 4 2 ...
##  $ TYPE_RESIDENCE     : Factor w/ 2 levels "PRINCIPALE","SECONDAIRE": 1 1 1 2 1 1 1 1 1 2 ...
##  $ TYPE_HABITATION    : Factor w/ 2 levels "APPARTEMENT",..: 1 2 1 2 2 2 2 2 2 1 ...
##  $ NB_PIECES          : int  1 NA 3 2 1 2 2 2 1 3 ...
##  $ SITUATION_JURIDIQUE: Factor w/ 2 levels "LOCATAIRE","PROPRIO": 2 2 1 1 1 1 1 1 2 1 ...
##  $ NIVEAU_JURIDIQUE   : Factor w/ 2 levels "JUR1","JUR2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ VALEUR_DES_BIENS   : num  3500 0 35000 9000 20000 9000 20000 9000 3500 0 ...
##  $ OBJETS_DE_VALEUR   : Factor w/ 2 levels "NIVEAU_1","NIVEAU_2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ZONIER             : Factor w/ 82 levels "A1","A10","A11",..: 48 3 39 75 82 71 70 81 14 12 ...
##  $ NBSIN_TYPE1_AN1    : int  0 0 0 0 0 0 1 0 0 0 ...
##  $ NBSIN_TYPE1_AN2    : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ NBSIN_TYPE1_AN3    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN1    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN2    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ NBSIN_TYPE2_AN3    : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ id                 : int  5 9 11 13 14 15 20 21 28 30 ...
##  $ ANNEE              : int  2017 2018 2017 2018 2017 2017 2018 2016 2018 2016 ...

Les données ont été prétraitées pour disposer d’une ligne par période de risque annuelle homogène.

Pour des raisons d’anonymisation, les identifiants contrats ont été supprimés de sorte qu’il n’est pas possible de suivre l’évolution du risque d’une période sur l’autre. Pour la même raison, il ne sera pas possible dans la construction des jeux de validation et de test de s’assurer que toute l’expérience d’un même assuré est bien intégralement dans un jeu d’entraînement, de validation ou de test.

La signification des colonnes présentes est la suivante :

X:

variable non nommée. Il s’agit d’un artefact du processus de création du fichier initial. On ne la prend pas en compte

EXPO :

exposition en année risque du contrat. Sa valeur précalculée pour chaque ligne est comprise entre 0 et 1

FORMULE :

variable catégorielle codant la formule du contrat comprenant 3 niveaux MEDIUM, ESSENTIEL et CONFORT

TYPE_RESIDENCE :

variable catégorielle codant le fait que le bien est une résidence PRINCIPALE, ou SECONDAIRE

TYPE_HABITATION :

variable catégorielle codant le fait que le bien est un APPARTEMENT, ou une MAISON

NB_PIECES :

variable numérique indiquant le nombre de pièces du logement

SITUATION_JURIDIQUE :

variable catégorielle indiquant si le souscripteur est prorpiétaire (PROPRIO) ou locataire (LOCATAIRE) du logement assuré

NIVEAU_JURIDIQUE :

variable catégorielle codant le niveau de couverture de la garantie juridique

VALEUR_DES_BIENS :

variable numérique reflétant la valeur couverte du contenu du logement

OBJETS_DE_VALEUR :

variable catégorielle codant le niveau de couverture des objets

ZONIER :

variable catégorielle codant la zone du bien assuré. Le code est constitué d’une lettre représentant une zone et d’un nombre représentant une partie de la zone

NBSIN_TYPE1_AN1 :

variable numérique indiquant le nombre de sinistre de type 1 l’année précédente

NBSIN_TYPE1_AN2 :

variable numérique indiquant le nombre de sinistre de type 1 il y a 2 ans

NBSIN_TYPE1_AN3 :

variable numérique indiquant le nombre de sinistre de type 1 il y a 3 ans

NBSIN_TYPE2_AN1 :

variable numérique indiquant le nombre de sinistre de type 2 l’année précédente

NBSIN_TYPE2_AN2 :

variable numérique indiquant le nombre de sinistre de type 2 il y a 2 ans

NBSIN_TYPE2_AN3 :

variable numérique indiquant le nombre de sinistre de type 2 il y a 3 ans

id :

identifiant du risque, clef de jointure

ANNEE :

variable catégorielle indiquant l’année d’observation du risque : 2016, 2017, 2018

3.2 sin_train

Le fichier expo_train.csv contient la liste des contrats en risques historiquement suivis :

sin_train = read.table(file = here('data', 'raw', 'sin_train.csv'), header=T, sep=';', dec=',', encoding = 'UTF-8', stringsAsFactors = T)

datatable(head(sin_train))

str(sin_train)

## 'data.frame':    2693 obs. of  4 variables:
##  $ id   : int  15 277 643 668 730 1816 1852 2061 2270 2322 ...
##  $ NB   : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ COUT : num  521.4 3000.2 26.3 462.4 640.9 ...
##  $ ANNEE: int  2017 2016 2018 2017 2018 2018 2016 2018 2018 2017 ...

id : identifiant du risque, clef de jointure avec le fichier expo_train
NB : Nombre de sinistres de la typologie modélisée (inconnue) sur la période d’observation

NB	nb_lignes
1	2,613.00
2	79.00
3	1.00

3.3 fichier combiné

On peut désormais combine les information d’exposition et de sinistres :

mrh <- expo_train %>%
      left_join(sin_train, by =c('id','ANNEE')) %>%
      replace_na(list('NB'=0, 'COUT'=0))

N’ayant pas accès au fichier de test utilisé lors du Hackathon, on dissocie ici le fichier mrh en un fichier train (70%), un fichier val (20%) et un fichier test (10%)

set.seed(42)

trainIndex<- createDataPartition(mrh$NB>=0, p=.7, list=FALSE, time=1)

mrh_train<-mrh[trainIndex,]
mrh_test<-mrh[ -trainIndex,]

valIndex<- createDataPartition(mrh_test$NB>=0, p=.66, list=FALSE, time=1)
mrh_val<-mrh_test[ valIndex,]
mrh_test<-mrh_test[ -valIndex,]

3.3.1 Taux de sinistres moyen :

tx <- sum(mrh_train$NB)/sum(mrh_train$EXPO)

print(paste0('taux de sinistre moyen annuel : ', round(tx * 100,2), '%'))

## [1] "taux de sinistre moyen annuel : 2.12%"

3.3.2 Taux de sinistre par variable :

3.3.2.1 FORMULE

FORMULE	Freq	nb_lignes	EXPO	pct_EXPO
ALL_INCLUDE	0.57%	4,695.00	3,355.67	3.67%
CONFORT	2.11%	61,372.00	51,858.85	56.79%
ESSENTIEL	2.22%	15,555.00	13,597.27	14.89%
MEDIUM	1.69%	27,334.00	22,501.75	24.64%

ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les niveaux CONFORT et ESSENTIEL semblent avoir un niveau de risque de fréquence équivalent

3.3.2.2 TYPE_RESIDENCE

TYPE_RESIDENCE	Freq	nb_lignes	EXPO	pct_EXPO
PRINCIPALE	2.15%	91,410.00	78,937.92	86.45%
SECONDAIRE	0.75%	17,546.00	12,375.63	13.55%

ENSEIGNEMENT : Le comportement observé est stable dans le temps. Les résidences SECONDAIRE ont systématiquement une fréquence moindre

3.3.2.3 TYPE_HABITATION

TYPE_HABITATION	Freq	nb_lignes	EXPO	pct_EXPO
APPARTEMENT	1.98%	49,749.00	43,093.59	47.19%
MAISON	1.95%	59,207.00	48,219.96	52.81%

ENSEIGNEMENT : Le comportement observé N’est PAS stable dans le temps. Le type d’habitation est à exclure.

3.3.2.4 NB_PIECES

NB_PIECES	Freq	nb_lignes	EXPO	pct_EXPO
0	0.12%	2,610.00	1,715.52	1.88%
1	1.61%	28,256.00	22,425.90	24.56%
2	1.83%	42,724.00	36,390.78	39.85%
3	2.48%	22,181.00	19,440.34	21.29%
4	2.95%	6,639.00	5,879.31	6.44%
NA	1.98%	6,546.00	5,461.69	5.98%

ENSEIGNEMENT : Le NB_PIECES présente des valeurs manquantes. La fréquence tend à augmenter avec le nombre de pièces.

3.3.2.5 SITUATION_JURIDIQUE

SITUATION_JURIDIQUE	Freq	nb_lignes	EXPO	pct_EXPO
LOCATAIRE	2.21%	60,339.00	53,225.98	58.29%
PROPRIO	1.63%	48,617.00	38,087.57	41.71%

ENSEIGNEMENT : Les propriétaires semblent avoir un fréquence moindre que les locataires de manière consistante, mais avec de la variabilité.

3.3.2.6 NIVEAU_JURIDIQUE

NIVEAU_JURIDIQUE	Freq	nb_lignes	EXPO	pct_EXPO
JUR1	1.97%	107,497.00	90,273.20	98.86%
JUR2	1.52%	1,459.00	1,040.35	1.14%

ENSEIGNEMENT : Le NIVEAU_JURIDIQUE N’est PAS stable dans le temps. A exclure.

3.3.2.7 VALEUR_DES_BIENS

VALEUR_DES_BIENS	Freq	nb_lignes	EXPO	pct_EXPO
0	0.25%	5,782.00	3,682.96	4.03%
3500	1.36%	36,543.00	28,888.32	31.64%
9000	1.95%	24,283.00	21,014.51	23.01%
20000	2.25%	23,639.00	20,768.83	22.74%
35000	2.86%	8,047.00	7,187.60	7.87%
50000	3.14%	5,470.00	4,991.12	5.47%
80000	3.09%	4,557.00	4,210.03	4.61%
100000	3.72%	635.00	570.18	0.62%

ENSEIGNEMENT : La fréquence augmente globalement avec la valeur des biens selon une relation non linéaire Par ailleurs au delà de 50k, l’expérience n’est pas stable. Les modalités supérieures à 80K devraient être groupées.

3.3.2.8 OBJETS_DE_VALEUR

OBJETS_DE_VALEUR	Freq	nb_lignes	EXPO	pct_EXPO
NIVEAU_1	1.84%	102,196.00	85,308.12	93.42%
NIVEAU_2	3.77%	6,760.00	6,005.44	6.58%

ENSEIGNEMENT : Le NIVEAU_2 présente une fréquence significativement plus élevée.

3.3.2.9 ZONIER

ZONIER	Freq	nb_lignes	EXPO	pct_EXPO
A1	1.75%	1,271.00	1,098.88	1.20%
A10	1.21%	768.00	659.36	0.72%
A11	0.62%	3,349.00	2,901.64	3.18%
A12	1.10%	2,182.00	1,896.53	2.08%
A13	0.97%	1,326.00	1,112.93	1.22%
A14	1.81%	596.00	497.59	0.54%
A2	1.14%	1,199.00	1,004.07	1.10%
A3	0.57%	814.00	705.81	0.77%
A4	0.00%	55.00	40.92	0.04%
A5	0.00%	11.00	10.24	0.01%
A6	0.72%	495.00	427.35	0.47%
A7	0.69%	170.00	145.53	0.16%
A8	2.29%	48.00	43.73	0.05%
A9	0.56%	2,885.00	2,488.88	2.73%
B1	0.41%	556.00	492.90	0.54%
B10	1.72%	2,492.00	2,140.03	2.34%
B11	0.00%	6.00	4.80	0.01%
B12	2.37%	249.00	210.74	0.23%
B13	6.23%	37.00	32.09	0.04%
B14	0.00%	83.00	69.97	0.08%
B16	0.76%	147.00	130.79	0.14%
B17	1.79%	2,442.00	2,071.38	2.27%
B18	0.00%	6.00	4.25	0.00%
B19	2.42%	672.00	572.48	0.63%
B2	1.09%	308.00	274.89	0.30%
B20	0.00%	6.00	5.58	0.01%
B21	0.00%	30.00	25.92	0.03%
B22	0.00%	41.00	36.37	0.04%
B23	0.78%	150.00	128.79	0.14%
B24	2.56%	91.00	78.21	0.09%
B25	1.33%	436.00	375.81	0.41%
B26	1.24%	1,031.00	884.64	0.97%
B27	0.84%	133.00	119.69	0.13%
B28	1.35%	173.00	148.10	0.16%
B29	1.83%	2,979.00	2,541.22	2.78%
B3	0.75%	1,394.00	1,210.68	1.33%
B30	0.96%	355.00	310.94	0.34%
B31	0.00%	56.00	51.74	0.06%
B32	1.81%	4,913.00	4,134.65	4.53%
B33	0.85%	134.00	117.90	0.13%
B34	1.23%	518.00	446.21	0.49%
B35	0.00%	3.00	2.28	0.00%
B36	1.33%	319.00	266.18	0.29%
B37	0.00%	108.00	94.32	0.10%
B38	5.38%	67.00	55.78	0.06%
B39	1.37%	1,804.00	1,542.57	1.69%
B4	0.73%	1,324.00	1,134.45	1.24%
B40	2.41%	4,045.00	3,427.84	3.75%
B41	2.79%	282.00	242.19	0.27%
B42	0.00%	22.00	20.55	0.02%
B43	1.72%	4,229.00	3,553.62	3.89%
B44	0.00%	5.00	4.09	0.00%
B45	0.83%	317.00	268.14	0.29%
B5	0.00%	24.00	20.47	0.02%
B6	1.30%	91.00	77.15	0.08%
B7	2.84%	243.00	211.00	0.23%
B8	0.00%	152.00	136.97	0.15%
B9	1.93%	437.00	374.18	0.41%
C1	2.68%	1,923.00	1,619.36	1.77%
C10	2.92%	638.00	515.30	0.56%
C11	3.19%	906.00	739.61	0.81%
C12	3.71%	1,700.00	1,335.15	1.46%
C13	2.09%	882.00	681.10	0.75%
C14	2.77%	962.00	769.95	0.84%
C15	3.70%	851.00	666.75	0.73%
C16	2.56%	894.00	707.78	0.78%
C17	2.18%	592.00	463.76	0.51%
C18	1.85%	526.00	429.41	0.47%
C19	3.29%	2,006.00	1,607.93	1.76%
C2	2.34%	4,216.00	3,597.31	3.94%
C20	3.06%	10,231.00	8,576.45	9.39%
C21	1.62%	3,365.00	2,789.43	3.05%
C22	1.50%	1,789.00	1,430.49	1.57%
C23	2.15%	4,418.00	3,642.23	3.99%
C24	1.91%	1,360.00	1,122.90	1.23%
C3	1.18%	2,127.00	1,817.00	1.99%
C4	3.82%	1,174.00	974.79	1.07%
C5	1.79%	4,772.00	3,952.35	4.33%
C6	1.59%	3,614.00	3,073.13	3.37%
C7	2.10%	1,442.00	1,195.34	1.31%
C8	1.65%	3,845.00	3,081.29	3.37%
C9	2.93%	6,644.00	5,438.74	5.96%

ENSEIGNEMENT : Le ZONIER présente une fréquence avec beaucoup de variabilité, mais il semble y avoir une croissance par région qu’on peut tester :

REGION	Freq	nb_lignes	EXPO	pct_EXPO
A	0.92%	15,169.00	13,033.47	14.27%
B	1.66%	32,910.00	28,052.53	30.72%
C	2.40%	60,877.00	50,227.55	55.01%

3.3.2.10 NBSIN_TYPE1_AN1

NBSIN_TYPE1_AN1	Freq	nb_lignes	EXPO	pct_EXPO
0	1.83%	101,793.00	84,856.05	92.93%
1	3.70%	6,581.00	5,945.96	6.51%
2	4.04%	530.00	469.28	0.51%
3	9.46%	52.00	42.26	0.05%

ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 l’année précédente. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.

3.3.2.11 NBSIN_TYPE1_AN2

NBSIN_TYPE1_AN2	Freq	nb_lignes	EXPO	pct_EXPO
1	1.85%	102,588.00	85,533.45	93.67%
2	3.50%	5,822.00	5,279.52	5.78%
3	4.92%	546.00	500.58	0.55%

ATTENTION : comme la table ci-dessous le démontre, cette variable présente un problème puisque toutes les lignes ont au moins un sinistre

Comme discuté dans le hackathon on ignore cette variable

3.3.2.12 NBSIN_TYPE1_AN3

NBSIN_TYPE1_AN3	Freq	nb_lignes	EXPO	pct_EXPO
0	1.86%	103,309.00	86,105.62	94.30%
1	3.42%	5,171.00	4,767.33	5.22%
2	6.13%	435.00	403.13	0.44%
3	5.34%	41.00	37.46	0.04%

ENSEIGNEMENT : La fréquence augmente globalement avec le nombre de sinistres de type 1 il y a 3 ans. Toutefois, il semble surtout q’il existe un effet d’une indicatrice a eu sinistre ou non.

3.3.2.13 NBSIN_TYPE2_AN1

NBSIN_TYPE2_AN1	Freq	nb_lignes	EXPO	pct_EXPO
0	1.95%	106,796.00	89,403.29	97.91%
1	2.61%	2,034.00	1,797.95	1.97%
2	2.93%	115.00	102.38	0.11%
3	0.00%	11.00	9.93	0.01%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.14 NBSIN_TYPE2_AN2

NBSIN_TYPE2_AN2	Freq	nb_lignes	EXPO	pct_EXPO
0	1.96%	95,382.00	79,861.28	87.46%
1	3.02%	1,452.00	1,307.09	1.43%
2	1.39%	77.00	71.91	0.08%
3	0.00%	3.00	3.00	0.00%
NA	1.86%	12,042.00	10,070.27	11.03%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.15 NBSIN_TYPE2_AN3

NBSIN_TYPE2_AN3	Freq	nb_lignes	EXPO	pct_EXPO
0	1.95%	107,445.00	89,931.57	98.49%
1	2.68%	1,414.00	1,292.12	1.42%
2	4.31%	93.00	85.86	0.09%
3	0.00%	4.00	4.00	0.00%

ENSEIGNEMENT : Il n’existe pas de lien stable sur la fréquence. A ne pas retenir.

3.3.2.16 ANNEE

ANNEE	Freq	nb_lignes	EXPO	pct_EXPO
2016	1.82%	37,947.00	31,635.56	34.64%
2017	2.04%	36,440.00	30,597.37	33.51%
2018	2.04%	34,569.00	29,080.62	31.85%

ENSEIGNEMENT : L’année 2016 semble atypique. Le niveau actuel semble plus proche de ceux de 2017 et 2018.

3.3.3 Nettoyage des données

preprocess<- function(dt) {
            res <- dt %>%
            mutate( 
              Y = NB/EXPO, #calcul de la fréquence pour XGB qui ne supporte pas les offset, mais accepte les poids
              NB_PIECES = ifelse(is.na(NB_PIECES), 2, NB_PIECES), #2 est le mode de NB_PIECES
              NBSIN_TYPE1_AN1 = ifelse(NBSIN_TYPE1_AN1==0,0,1), #On regroupe les modalités non stables dans le temps
              NBSIN_TYPE1_AN3 = ifelse(NBSIN_TYPE1_AN3==0,0,1), #On regroupe les modalités non stables dans le temps
              REGION =  gsub('[0-9]', '', ZONIER), #création de région
              VALEUR_DES_BIENS = pmin(VALEUR_DES_BIENS, 50000), # On regroupe les valeurs des biens supérieurs à 80k
              OFFSET = log(EXPO) # création d'une variable d'offset qui servira pour le gbm
            ) %>%
            mutate( REGION = as.factor(REGION),
                    VALEUR_DES_BIENS = as.factor(VALEUR_DES_BIENS)
            ) %>%
            select(-c('X', 'TYPE_HABITATION', 'NIVEAU_JURIDIQUE', 'NBSIN_TYPE1_AN2', 'NBSIN_TYPE2_AN1', 'NBSIN_TYPE2_AN2', 'NBSIN_TYPE2_AN3', 'ZONIER'  )) %>% # On supprime les variables sans lien stable dans le temps
            relocate('id', 'EXPO', 'FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE', 'OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES',  'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE', 'NB', 'Y', 'COUT' )

            return(res)
}

df_train <- preprocess(mrh_train)
df_val <- preprocess(mrh_val)
df_test <- preprocess(mrh_test)

3.3.4 Calcul du V de Cramer

On mesude l’association des variables entre elles

X_train <- df_train %>% select(-c('id', 'EXPO', 'NB', 'Y', 'COUT'))


# Initialize empty matrix to store coefficients
empty_m <- matrix(ncol = length(X_train),
            nrow = length(X_train),
            dimnames = list(names(X_train), 
                            names(X_train)))
# Function that accepts matrix for coefficients and data and returns a correlation matrix
calculate_cramer <- function(m, df) {
 for (r in seq(nrow(m))){
   for (c in seq(ncol(m))){
     m[[r, c]] <- lsr::cramersV(X_train[[r]], X_train[[c]])
   }
 }
    return(m)
}

cor_matrix <- calculate_cramer(empty_m ,X_train)

corrplot(cor_matrix)

4 Modélisation de la fréquence

Dans cette partie, on essaye de modéliser la fréquence des sinistres.

summary(df_test$NB)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01745 0.00000 2.00000

4.1 Modèle de référence

On commence par une modélisation de référence réduite à un modèle GLM POisson avec un unique intercept.

reg0 = glm(NB~1+offset(OFFSET),family=poisson,data=df_train)

summary(reg0)

## 
## Call:
## glm(formula = NB ~ 1 + offset(OFFSET), family = poisson, data = df_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.2058  -0.2058  -0.2058  -0.1774   4.8750  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -3.85523    0.02274  -169.5   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 15507  on 108955  degrees of freedom
## Residual deviance: 15507  on 108955  degrees of freedom
## AIC: 19294
## 
## Number of Fisher Scoring iterations: 6

4.2 Modèle GLM plus élaboré

Il s’agit d’un modèle GLM Poisson classique tirant partie des enseignements de la partie précédente, sans interaction. Les variables sont traitées comme catégorielles pour permettre les effets non linéaires.

fit_glm <- glm(NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE), 
              data=df_train, offset=OFFSET, family=quasipoisson())

summary(fit_glm)

## 
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + SITUATION_JURIDIQUE + 
##     OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + 
##     REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + 
##     as.factor(ANNEE), family = quasipoisson(), data = df_train, 
##     offset = OFFSET)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4898  -0.2195  -0.1769  -0.1220   4.3196  
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -8.35103    0.78732 -10.607  < 2e-16 ***
## FORMULECONFORT                    0.89291    0.21816   4.093 4.26e-05 ***
## FORMULEESSENTIEL                  0.85630    0.22564   3.795 0.000148 ***
## FORMULEMEDIUM                     0.85559    0.21928   3.902 9.55e-05 ***
## TYPE_RESIDENCESECONDAIRE         -0.61485    0.11032  -5.573 2.51e-08 ***
## SITUATION_JURIDIQUEPROPRIO       -0.05746    0.06009  -0.956 0.338984    
## OBJETS_DE_VALEURNIVEAU_2          0.26778    0.08224   3.256 0.001129 ** 
## as.factor(VALEUR_DES_BIENS)3500   0.62720    0.34944   1.795 0.072676 .  
## as.factor(VALEUR_DES_BIENS)9000   0.77173    0.35190   2.193 0.028306 *  
## as.factor(VALEUR_DES_BIENS)20000  0.83537    0.35252   2.370 0.017803 *  
## as.factor(VALEUR_DES_BIENS)35000  0.97892    0.35803   2.734 0.006254 ** 
## as.factor(VALEUR_DES_BIENS)50000  0.93078    0.35814   2.599 0.009352 ** 
## as.factor(NB_PIECES)1             2.04164    0.82753   2.467 0.013621 *  
## as.factor(NB_PIECES)2             1.98341    0.82634   2.400 0.016387 *  
## as.factor(NB_PIECES)3             2.23041    0.82810   2.693 0.007074 ** 
## as.factor(NB_PIECES)4             2.31557    0.83137   2.785 0.005350 ** 
## REGIONB                           0.43688    0.10543   4.144 3.42e-05 ***
## REGIONC                           0.94976    0.09974   9.522  < 2e-16 ***
## as.factor(NBSIN_TYPE1_AN1)1       0.43828    0.07264   6.034 1.61e-09 ***
## as.factor(NBSIN_TYPE1_AN3)1       0.33205    0.08158   4.070 4.70e-05 ***
## as.factor(ANNEE)2017              0.10332    0.05889   1.754 0.079351 .  
## as.factor(ANNEE)2018              0.10273    0.05960   1.724 0.084798 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.119556)
## 
##     Null deviance: 15507  on 108955  degrees of freedom
## Residual deviance: 14906  on 108934  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 9

Le modèle n’explique qu’une faible partie de la variabilité.

Par ailleurs la variable SITUATION_JURIDIQUE (qui n’a que 2 niveaux) n’a finalement pas d’impact sur le risque et peut donc être exclue du modèle :

fit_glm <- glm(NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR + as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + as.factor(ANNEE), 
              data=df_train, offset=OFFSET, family=quasipoisson())

summary(fit_glm)

## 
## Call:
## glm(formula = NB ~ FORMULE + TYPE_RESIDENCE + OBJETS_DE_VALEUR + 
##     as.factor(VALEUR_DES_BIENS) + as.factor(NB_PIECES) + REGION + 
##     as.factor(NBSIN_TYPE1_AN1) + as.factor(NBSIN_TYPE1_AN3) + 
##     as.factor(ANNEE), family = quasipoisson(), data = df_train, 
##     offset = OFFSET)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.4903  -0.2202  -0.1776  -0.1224   4.3204  
## 
## Coefficients:
##                                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                      -8.41067    0.78465 -10.719  < 2e-16 ***
## FORMULECONFORT                    0.90472    0.21774   4.155 3.25e-05 ***
## FORMULEESSENTIEL                  0.87116    0.22504   3.871 0.000108 ***
## FORMULEMEDIUM                     0.86290    0.21908   3.939 8.20e-05 ***
## TYPE_RESIDENCESECONDAIRE         -0.60328    0.10960  -5.504 3.71e-08 ***
## OBJETS_DE_VALEURNIVEAU_2          0.26617    0.08219   3.238 0.001203 ** 
## as.factor(VALEUR_DES_BIENS)3500   0.61153    0.34890   1.753 0.079655 .  
## as.factor(VALEUR_DES_BIENS)9000   0.76481    0.35165   2.175 0.029640 *  
## as.factor(VALEUR_DES_BIENS)20000  0.83693    0.35233   2.375 0.017530 *  
## as.factor(VALEUR_DES_BIENS)35000  0.98370    0.35781   2.749 0.005974 ** 
## as.factor(VALEUR_DES_BIENS)50000  0.93646    0.35790   2.616 0.008885 ** 
## as.factor(NB_PIECES)1             2.06378    0.82697   2.496 0.012577 *  
## as.factor(NB_PIECES)2             2.01424    0.82547   2.440 0.014684 *  
## as.factor(NB_PIECES)3             2.27101    0.82677   2.747 0.006018 ** 
## as.factor(NB_PIECES)4             2.35781    0.82997   2.841 0.004500 ** 
## REGIONB                           0.43567    0.10539   4.134 3.57e-05 ***
## REGIONC                           0.94502    0.09959   9.489  < 2e-16 ***
## as.factor(NBSIN_TYPE1_AN1)1       0.44183    0.07254   6.091 1.13e-09 ***
## as.factor(NBSIN_TYPE1_AN3)1       0.33622    0.08146   4.128 3.67e-05 ***
## as.factor(ANNEE)2017              0.10379    0.05887   1.763 0.077911 .  
## as.factor(ANNEE)2018              0.10351    0.05958   1.737 0.082354 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for quasipoisson family taken to be 1.119014)
## 
##     Null deviance: 15507  on 108955  degrees of freedom
## Residual deviance: 14907  on 108935  degrees of freedom
## AIC: NA
## 
## Number of Fisher Scoring iterations: 9

fit_glm_Diag <- data.frame(df_test,
                     fit = predict(fit_glm, newdata = df_test, type="response")
                     ) %>%
                     mutate(pearson = (NB-fit)/sqrt(fit))

dat <- fit_glm_Diag %>%
       select( fit, pearson) %>%
       mutate(quantile = ntile(fit, 16000)) %>%
       group_by(quantile) %>%
       summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))


ggplot(dat, aes(x = mean_fit, y = mean_pearson)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 0, size = 1, color="red") + 
  ggtitle("Rédidus de Pearson aggrégés - test - GLM")

round(sum(fit_glm_Diag$fit)/sum(df_test$EXPO),4)

## [1] 0.0214

round(sum(df_test$NB)/sum(df_test$EXPO),4)

## [1] 0.0209

fit_glm_Diag %>% 
  mutate(residual = NB-fit) %>%
  summarize(rmse  = sqrt(mean(residual^2)))

##        rmse
## 1 0.1329575

summary(df_test$NB)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.00000 0.00000 0.01745 0.00000 2.00000

summary(fit_glm_Diag$fit)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.440e-06 8.047e-03 1.632e-02 1.784e-02 2.448e-02 1.162e-01

clrpal4 = wes_palette("Darjeeling1")[c(1,2,3,5)]


hist(fit_glm_Diag$fit, breaks = 50, xlim=c(0,.08),col=clrpal4[1],border="white",
     main="",xlab="Premium (GLM)",ylab="",ylim=c(0,2000))

4.3 Modèle de Gradient Boosting

Pour la mise en place d’un modèle GBM, nous utiliserons le package h2o :

h2o.init(nthreads = -1)

h2o.no_progress()

Ce package permet de définir une fonction de perte correspondant à la déviance de Poisson avec sa fonction de lien canonique, c’est à dire le ln dans notre cas. Par ailleurs, la méthode permet de définir une colonne qui peut servir d’offset. D’après la documentation, cet offset est appliqué dans l’espace linéarisé :

For other distributions, the offset corrections are applied in the linearized space before applying the inverse link function to get the actual response values.

Il nous faut donc passer en valeur le $ln(EXPO)$ que nous avons renseigné dans la variable OFFSET :

fit_gbm = h2o.gbm(y = 'NB', x = names(df_train)[3:12],
                       distribution = "poisson",
                       offset_column = "OFFSET",
                       stopping_metric = "deviance",
                       training_frame = as.h2o(df_train),
                       validation_frame = as.h2o(df_val),
                       ntrees = 60,
                       nfolds = 5,
                       seed = 1)

## Depth 10 is usually plenty of depth for most datasets, but you never know
hyper_params = list( max_depth = seq(1,15,2) )
#hyper_params = list( max_depth = c(4,6,8,12,16,20) ) ##faster for larger datasets


grid <- h2o.grid(
  ## hyper parameters
  hyper_params = hyper_params,
  ## full Cartesian hyper-parameter search
  search_criteria = list(strategy = "Cartesian"),
  ## which algorithm to run
  algorithm="gbm",
  ## identifier for the grid, to later retrieve it
  grid_id="Grid_gbm_offset",
  ## standard model parameters
  distribution = "poisson",
  x = names(df_train)[3:12],
  y = 'NB',
  offset_column = "OFFSET",
  training_frame = as.h2o(df_train),
  validation_frame = as.h2o(df_val),

  ## more trees is better if the learning rate is small enough
  ## here, use "more than enough" trees - we have early stopping
  ntrees = 10000,
  ## smaller learning rate is better
  ## since we have learning_rate_annealing, we can afford to start with a bigger learning rate
  learn_rate = 0.05,
  ## learning rate annealing: learning_rate shrinks by 1% after every tree
  ## (use 1.00 to disable, but then lower the learning_rate)
  learn_rate_annealing = 0.99,
  ## sample 80% of rows per tree
  sample_rate = 0.8,
  ## sample 80% of columns per split
  col_sample_rate = 0.8,
  ## fix a random number generator seed for reproducibility
  seed = 1234,
  ## early stopping once the validation AUC doesn't improve by at least 0.01% for 5 consecutive scoring events
  stopping_rounds = 5,
  stopping_tolerance = 1e-4,
  stopping_metric = "deviance",
  ## score every 10 trees to make early stopping reproducible (it depends on the scoring interval)
  score_tree_interval = 10
)
## by default, display the grid search results sorted by increasing logloss (since this is a classification task)
grid

## H2O Grid Details
## ================
## 
## Grid ID: Grid_gbm_offset 
## Used hyper parameters: 
##   -  max_depth 
## Number of models: 8 
## Number of failed models: 0 
## 
## Hyper-Parameter Search Summary: ordered by increasing residual_deviance
##   max_depth               model_ids   residual_deviance
## 1         3 Grid_gbm_offset_model_2 0.17448113788335165
## 2         1 Grid_gbm_offset_model_1 0.17561652746681827
## 3         5 Grid_gbm_offset_model_3  0.1798611385902214
## 4         7 Grid_gbm_offset_model_4 0.19916164139992712
## 5         9 Grid_gbm_offset_model_5   0.246858307237967
## 6        11 Grid_gbm_offset_model_6  0.3093800261468492
## 7        13 Grid_gbm_offset_model_7  0.3404180725614419
## 8        15 Grid_gbm_offset_model_8  0.3483873140442045

fit_gbm <- h2o.getModel(grid@model_ids[[1]])

print(h2o.performance(fit_gbm, newdata = as.h2o(df_val)))

## H2ORegressionMetrics: gbm
## 
## MSE:  0.01898311
## RMSE:  0.1377792
## MAE:  0.03495651
## RMSLE:  0.0935029
## Mean Residual Deviance :  0.1744811

plot(fit_gbm)

h2o.performance(fit_gbm, newdata = as.h2o(df_test))

## H2ORegressionMetrics: gbm
## 
## MSE:  0.0176801
## RMSE:  0.1329665
## MAE:  0.03420818
## RMSLE:  0.0912325
## Mean Residual Deviance :  0.1701883

fit_gbm_Diag <- data.frame(df_test,
                     fit = as.vector(h2o.predict(object = fit_gbm, newdata=as.h2o(df_test)))
                     ) %>%
                     mutate(pearson = (NB-fit)/sqrt(fit))


dat_gbm <- fit_gbm_Diag %>%
       select( fit, pearson) %>%
       mutate(quantile = ntile(fit, 16000)) %>%
       group_by(quantile) %>%
       summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))


ggplot(dat_gbm, aes(x = mean_fit, y = mean_pearson)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 0, size = 1, color="red") + 
  ggtitle("Rédidus de Pearson aggrégés - test - GBM")

round(sum(fit_gbm_Diag$fit)/sum(df_test$EXPO),4)

## [1] 0.0211

round(sum(df_test$NB)/sum(df_test$EXPO),4)

## [1] 0.0209

fit_gbm_Diag %>% 
  mutate(residual = NB-fit) %>%
  summarize(rmse  = sqrt(mean(residual^2)))

##        rmse
## 1 0.1329665

summary(fit_gbm_Diag$fit)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 7.950e-06 7.993e-03 1.675e-02 1.759e-02 2.453e-02 1.084e-01

hist(fit_gbm_Diag$fit, breaks = 100, xlim=c(0,.08),col=clrpal4[2],border="white",
     main="",xlab="Premium (GBM)",ylab="",ylim=c(0,4000))

4.4 Modèle de XGBoost

x <- c('FORMULE', 'TYPE_RESIDENCE', 'SITUATION_JURIDIQUE','OBJETS_DE_VALEUR', 'VALEUR_DES_BIENS', 'NB_PIECES', 'REGION', 'NBSIN_TYPE1_AN1', 'NBSIN_TYPE1_AN3', 'ANNEE')
y <- 'Y'
w <- 'EXPO'

# Input maker
prep_xgb <- function(dat, x) {
  data.matrix(dat[, x, drop = FALSE])
}
# Data interface to XGBoost
dtrain <- xgb.DMatrix(
  prep_xgb(df_train, x), 
  label = df_train[[y]], 
  weight = df_train[[w]]
)
# Parameters chosen by 5-fold grouped CV
params_freq <- list(
  learning_rate = 0.2,
  max_depth = 5,
  alpha = 3,
  lambda = 0.5,
  max_delta_step = 2,
  min_split_loss = 0,
  colsample_bytree = 1,
  subsample = 0.9
)



# Fit
set.seed(1)
fit_xgb <- xgb.train(
  params_freq, 
  data = dtrain,
  nrounds = 580,
  objective = "count:poisson",
  watchlist = list(train = dtrain),
  print_every_n = 100
)

## [1]  train-poisson-nloglik:0.504222 
## [101]    train-poisson-nloglik:0.114902 
## [201]    train-poisson-nloglik:0.100886 
## [301]    train-poisson-nloglik:0.100259 
## [401]    train-poisson-nloglik:0.099847 
## [501]    train-poisson-nloglik:0.099479 
## [580]    train-poisson-nloglik:0.099241

# Save and load model
xgb.save(fit_xgb, "xgb.model")

## [1] TRUE

fit_xgb <- xgb.load("xgb.model")

fit_xgb_Diag <- data.frame(df_test,
                     fit = predict(fit_xgb, prep_xgb(df_test, x))
                     ) %>%
                     mutate(pearson = (NB-fit)/sqrt(fit))


dat_xgb <- fit_xgb_Diag %>%
       select( fit, pearson) %>%
       mutate(quantile = ntile(fit, 16000)) %>%
       group_by(quantile) %>%
       summarize(mean_fit = mean(fit, na.rm = TRUE), mean_pearson = mean(pearson, na.rm = TRUE))


ggplot(dat_xgb, aes(x = mean_fit, y = mean_pearson)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 0, size = 1, color="red") + 
  ggtitle("Rédidus de Pearson aggrégés - test - XGB")

round(sum(fit_xgb_Diag$fit)/sum(df_test$EXPO),4)

## [1] 0.0249

round(sum(df_test$NB)/sum(df_test$EXPO),4)

## [1] 0.0209

fit_xgb_Diag %>% 
  mutate(residual = NB-fit) %>%
  summarize(rmse  = sqrt(mean(residual^2)))

##        rmse
## 1 0.1330936

hist(fit_xgb_Diag$fit, breaks = 200, xlim=c(0,.08),col=clrpal4[3],border="white",
     main="",xlab="Premium (GBM)",ylab="",ylim=c(0,4000))

5 Model Explanations

The models are ready, so let’s shed light into them.

5.1 Setting up explainers

We start by setting up the explainers. These are basically objects that know how to create predictions.

Crucial: the prediction function needs to work for subsets of datasets, not just for the original model data set.

fl_glm <- flashlight(
  model = fit_glm, label = "GLM", 
  predict_function = function(fit, X) predict(fit, X, type = "response")
)

fl_xgb <- flashlight(
  model = fit_xgb, label = "XGBoost", 
  predict_function = function(fit, X) predict(fit, prep_xgb(X, x))
)

fl_gbm <- flashlight(
  model = fit_gbm, label = "H20 GBM", 
  predict_function = function(fit, X) as.vector(h2o.predict(object = fit_gbm, newdata=as.h2o(X)))
)


# Combine them and add common elements like reference data
metrics <- list(`Average deviance` = deviance_poisson, 
                `Relative deviance reduction` = r_squared_poisson)
fls <- multiflashlight(list(fl_glm, fl_xgb, fl_gbm), data = df_test, 
                       y = y, w = w, metrics = metrics)
# Version on canonical scale
fls_log <- multiflashlight(fls, linkinv = log)

5.2 Performance

We start to interpret the models by looking at model performance using deviance related metrics.

fillc <- "#E69F00"
perf <- light_performance(fls)
perf

## 
## I am an object of class light_performance_multi 
## 
## data.frames:
## 
##  data 
## # A tibble: 6 x 3
##   label   metric                       value
##   <fct>   <fct>                        <dbl>
## 1 GLM     Average deviance            0.166 
## 2 GLM     Relative deviance reduction 0.0107
## 3 XGBoost Average deviance            0.163 
## 4 XGBoost Relative deviance reduction 0.0280
## 5 H20 GBM Average deviance            0.166 
## 6 H20 GBM Relative deviance reduction 0.0103

plot(perf, geom = "point") +
  labs(x = element_blank(), y = element_blank())

5.3 Importance

Next, we consider permutation variable importance. By default flashlight uses the first performance metric specified above.

imp <- light_importance(fls, v = x)
plot(imp, fill = fillc, color = "black")

5.3.1 Partial dependence curves

Taking the average of many ICE curves produces the famous partial dependence plot. We consider such plots for four predictors.

plot(light_profile(fls, v = "VALEUR_DES_BIENS"))

plot(light_profile(fls, v = "NB_PIECES"))

plot(light_profile(fls, v = "REGION"))

5.4 Classic diagnostic plots

Classic predicted/residual/response versus covariable plots are worth a look.

# Average predicted versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "predicted"))

# Average residual versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "residual")) +
  geom_hline(yintercept = 0)

# Average response versus covariable
plot(light_profile(fls, v = "NB_PIECES", type = "response"))

5.4.1 Multiple aspects combined

We often get a good picture of the effect of a covariable by combining partial dependence with classic plots.

eff_NB_PIECES <- light_effects(fls, v = "NB_PIECES", counts_weighted = TRUE)
p <- plot(eff_NB_PIECES, show_points = FALSE)
plot_counts(p, eff_NB_PIECES, alpha = 0.3)

interact_rel <- light_interaction(
  fls_log, 
  v = most_important(imp, 4), 
  take_sqrt = FALSE,
  pairwise = TRUE, 
  use_linkinv = TRUE,
  seed = 61
)
plot(interact_rel, color = "black", fill = fillc, rotate_x = TRUE)

5.4.2 On absolute scale

interact_abs <- light_interaction(
  fls_log, 
  v = most_important(imp, 4), 
  normalize = FALSE,
  pairwise = TRUE, 
  use_linkinv = TRUE,
  seed = 61
)
plot(interact_abs, color = "black", fill = fillc, rotate_x = TRUE)

5.4.3 Visualization

Let’s illustrate a strong and a weak interaction by conditional partial dependence plots.

# Strong interaction
# p <- plot(eff_DrivAge, show_points = FALSE)
# plot_counts(p, eff_DrivAge, alpha = 0.3)


pdp_VALEUR_DES_BIENS_REGION <- light_profile(fls_log, v = "VALEUR_DES_BIENS", by = "REGION", 
                                  pd_seed = 50, data = df_val, counts_weighted=TRUE)
p <- plot(pdp_VALEUR_DES_BIENS_REGION, show_points = FALSE)
pdp_count <- light_effects(fls, v = "REGION")
plot_counts(p, pdp_count, alpha = 0.3)

pdp_VALEURDESBIENS_TYPE_RESIDENCE <- light_profile(fls_log, v = "VALEUR_DES_BIENS", by = "TYPE_RESIDENCE", 
                                  pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_VALEURDESBIENS_TYPE_RESIDENCE)

# Weak interaction
pdp_TYPE_RESIDENCE_NB_PIECES <- light_profile(fls_log, v = "TYPE_RESIDENCE", 
                                 by = "NB_PIECES", pd_seed = 50, data = df_val, counts_weighted=TRUE)
plot(pdp_TYPE_RESIDENCE_NB_PIECES)

5.5 Primes

Comme on l’a vu à la section précédente, le jeu de données est assez simple et offre peu d’opportunités de prise en compte d’effets non linéaires. Néanmoins même dans ces conditions on peut observer que si le modèles GLM et GBM sont assez proches de la référence observée sur le portefeuille, on observe que le XGBoost surestime dans ce cas de presque 19% la charge au global :

dat = df_test

# print (paste0('fréquence de référence sur le jeu de données : ', sum(dat$NB)/sum(dat$EXPO)))
# print (paste0('fréquence estimée par le GLM sur le jeu de données : ', sum(fl_glm$predict_function(fit_glm, dat))/sum(dat$EXPO)))
# print (paste0('fréquence estimée par le GBM sur le jeu de données : ', sum(fl_gbm$predict_function(fit_gbm, dat))/sum(dat$EXPO)))
# print (paste0('fréquence estimée par le XGB sur le jeu de données : ', sum(fl_xgb$predict_function(fit_xgb, dat))/sum(dat$EXPO)))

ref <- data.frame(modeles=c('référence', 'GLM', 'GBM', 'XGB'))
ref$primesmoy = c(sum(dat$NB)/sum(dat$EXPO),
                  sum(fl_glm$predict_function(fit_glm, dat))/sum(dat$EXPO),
                  sum(fl_gbm$predict_function(fit_gbm, dat))/sum(dat$EXPO),
                  sum(fl_xgb$predict_function(fit_xgb, dat))/sum(dat$EXPO)
)
ref <- ref%>%
       mutate(PctECART = primesmoy/ (sum(dat$NB)/sum(dat$EXPO)) -1 ) %>%
       mutate(PctECART = sprintf("%.1f %%", 100*PctECART) )
ref

##     modeles  primesmoy PctECART
## 1 référence 0.02090040    0.0 %
## 2       GLM 0.02137483    2.3 %
## 3       GBM 0.02107068    0.8 %
## 4       XGB 0.02485922   18.9 %

dat = df_val
valid_prime_glm = fl_glm$predict_function(fit_glm, dat)
valid_prime_gbm = fl_gbm$predict_function(fit_gbm, dat)
valid_prime_xgb = fl_xgb$predict_function(fit_xgb, dat)


par(mfrow=c(1,3))
hist(valid_prime_glm, breaks = c((0:(201*max(valid_prime_glm)))/200, 201/200*max(valid_prime_glm)),xlim=c(0,.2),col=clrpal4[1],border="white",
     main="",xlab="Primes estimées (GLM)",ylab="",ylim=c(0,37000))
hist(valid_prime_gbm, breaks = c((0:(201*max(valid_prime_gbm)))/200, 201/200*max(valid_prime_gbm)),xlim=c(0,.2),col=clrpal4[2],border="white",
     main="",xlab="Primes estimées (GBM)",ylab="",ylim=c(0,37000))
hist(valid_prime_xgb, breaks = c((0:(201*max(valid_prime_xgb)))/200, 201/200*max(valid_prime_xgb)),xlim=c(0,.2),col=clrpal4[3],border="white",
     main="",xlab="Primes estimées (XGB)",ylab="",ylim=c(0,37000))

DSA 2021 - TD R

Fabien FAIVRE

18 juillet 2021

1 Présentation du projet

1.1 Les bonnes propriétés du GLM dans le cadre d’une modélisation de la fréquence de sinistres

1.2 L’apport et les limites des modèles ML en tarification

1.2.1 Ajustement #1 : définition de la fonction de perte

1.2.2 Limites de la modification de la fonction de perte

1.3 Pistes de réconciliation des approches GLM et ML pour la tarification

1.4 Objectif du projet

2 Analyse

2.1 Mise en place de l’environnement

2.2 code

3 Présentation du jeu de données

3.1 expo_train

3.2 sin_train

3.3 fichier combiné

3.3.1 Taux de sinistres moyen :

3.3.2 Taux de sinistre par variable :

3.3.2.1 FORMULE

3.3.2.2 TYPE_RESIDENCE

3.3.2.3 TYPE_HABITATION

3.3.2.4 NB_PIECES

3.3.2.5 SITUATION_JURIDIQUE

3.3.2.6 NIVEAU_JURIDIQUE

3.3.2.7 VALEUR_DES_BIENS

3.3.2.8 OBJETS_DE_VALEUR

3.3.2.9 ZONIER

3.3.2.10 NBSIN_TYPE1_AN1

3.3.2.11 NBSIN_TYPE1_AN2

3.3.2.12 NBSIN_TYPE1_AN3

3.3.2.13 NBSIN_TYPE2_AN1

3.3.2.14 NBSIN_TYPE2_AN2

3.3.2.15 NBSIN_TYPE2_AN3

3.3.2.16 ANNEE

3.3.3 Nettoyage des données

3.3.4 Calcul du V de Cramer

4 Modélisation de la fréquence

4.1 Modèle de référence

4.2 Modèle GLM plus élaboré

4.3 Modèle de Gradient Boosting

4.4 Modèle de XGBoost

5 Model Explanations

5.1 Setting up explainers

5.2 Performance

5.3 Importance

5.3.1 Partial dependence curves

5.4 Classic diagnostic plots

5.4.1 Multiple aspects combined

5.4.2 On absolute scale

5.4.3 Visualization

5.5 Primes